%pylab inline
import scipy
import scipy.signal
import scipy.fftpack
import numpy as np

Fs=962 # Sample rate, Hz

Populating the interactive namespace from numpy and matplotlib

RUN_NAME="candle"

if False:
    import serial

    s = serial.Serial("/dev/tty.usbserial-A4006Dx4", 9600)

    N_SAMPLES=Fs*60


    data = [ ord(s.read()) for _ in range(N_SAMPLES) ]

    s.close()

    np.save("%s.npy", data)
else:
    data = np.load("%s.npy" % RUN_NAME)
    
figure(figsize=(20,5))
plot(data)
len(data)

57720

def show_fft(Fs, y):
    # sample spacing
    T = 1.0 / Fs
    # Number of samplepoints
    N = len(y)

    w = scipy.signal.blackman(N)
    ywf = scipy.fftpack.fft(y*w)

    xf = np.linspace(0.0, 1.0/(2.0*T), N/2)

#    axvline(60, color="#ffcccc")
#    axvline(120, color="#ccccff")


    semilogy(xf[1:N/2], 2.0/N * np.abs(ywf[1:N/2]), '-r')
    grid()
  

def show_fft_and_plot(Fs, y):
    figure(figsize=(20,5))

    subplot(1,2,1)
    show_fft(Fs,y)

    subplot(1,2,2)
    plot(y)

t0 = 5700
dt = 5000

show_fft_and_plot(Fs, data[t0:t0+dt])

lpf = scipy.signal.remez(150,[ 0.0,0.09,0.15,0.5], [1,0])
yp = scipy.signal.lfilter(lpf, [1], y)



show_fft_and_plot(Fs, yp)

figure(figsize=(20,5))
hist(yp, bins=100)
ylim(0,400)
True

True

hpf = scipy.signal.remez(147,[ 0.0,0.2, 0.25,0.5], [0,1])
yp = scipy.signal.lfilter(hpf, [1], y)

show_fft_and_plot(Fs, yp)

figure(figsize=(20,5))
hist(yp, bins=100)
ylim(0,200)
True

figure(figsize=(20,5))
show_fft(Fs, yp)
ylim(0.001,2)

(0.001, 2)

def gen_f(Fs, f, n):
    return np.sin(2*math.pi*f/Fs*np.linspace(0.0,n,n))


def downmix(Fs, f, y):
    return y*gen_f(Fs, f, len(y))

figure(figsize=(20,5))
#plot()

lpf_am_demod = scipy.signal.remez(50,[0,0.03,0.05,0.5], [1,0])

yd = downmix(Fs, 140, diff(y))
y_demod = scipy.signal.lfilter(lpf_am_demod, [1], abs(yd) )
#y_demod = scipy.signal.lfilter(lpf_am_demod, [1], abs(downmix(Fs, 280, yp)))  + \
#          scipy.signal.lfilter(lpf_am_demod, [1], abs(downmix(Fs, 140, yp)))  + \
#          scipy.signal.lfilter(lpf_am_demod, [1], abs(downmix(Fs, 420, yp)) )

plot(diff(y), alpha=0.2)
plot(y_demod[50:]*2, 'k')


ylim(0,300)

(0, 300)